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Abstract 


Nonlinear time-dependent calculations have been carried out in order to study the 
evolution of the thermal instability for vertically integrated, non self-gravitating 
models of optically thick, transonic, slim accretion discs around black holes. In these 
calculations we investigated only the original version of the slim disc model with low 
viscosity {a = 0.001) and for a stellar mass (10 Mq) black hole. This original version 
of the model does not yet include several important non-local effects (viscous forces 
in the radial equation of motion, diffusion-type formulation for the viscosity in the 
angular momentum equation, viscous dissipation rate associated with the stress in 
the azimuthal direction and radiative losses in the radial direction in the energy 
balance equation). It is clear, therefore, that this treatment is greatly simplihed but 
our strategy is to consider this as a standard reference against which to compare 
results from forthcoming studies in which the additional effects will be added one 
by one thus giving a systematic way of understanding the contribution from each of 
them. 

We have considered an already-formed disc and studied its stability against small 
axisymmetric perturbations. Those models which were stable according to local 
analysis, remain stable and stationary to a good approximation as also do models 
for which local analysis predicted an unstable region with radial dimension smaller 
than the shortest wavelength of the unstable modes. In terms of luminosity, all 
models with luminosity less than or equal to 0.08L^ are stable. For models with 
higher luminosity than this but which are still sub-Eddington, a shock-like structure 
forms near to the sonic point, probably leading to subsequent disruption of the 
disc. The model with Eddington luminosity evolves in a violent way, with the 
shock-like structure being formed already within the hrst second of evolution. From 
an examination of phase space trajectories, our preliminary conclusion is that the 
stabilizing effect of advection is not strong enough in these models to allow for limit 
cycle behaviour to occur. However, in order to make a dehnitive statement on this, 
it would be necessary to implement special numerical techniques for treatment of 
the shock-like structure, which becomes very extreme, and this lies beyond the scope 
of the present paper. 

Subject headings: accretion: accretion discs, instabilities: thermal 
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1 Introduction 


Accretion discs around black holes are thought to be present in galactic black hole 
candidates and active galactic nuclei. One of the characteristic properties of these 
objects is their variability. Apart from fundamental questions concerning the very 
existence of fluid disc-like conhgurations and their lifetimes, the most immediate 
purpose of the study of accretion disc stability is to discover the origin of the in¬ 
cipient instabilities which may be responsible for the observed variability. However, 
before applying the results of stability studies to realistic objects, one must elimi¬ 
nate the possibility that the behaviour of the model is merely a consequence of the 
approximations made in its construction (Piran, 1978; Pringle, 1981). 

Equilibrium models of geometrically thin accretion discs were found to be subject 
to two types of instability - thermal (Pringle, Rees, and Pacholczyk, 1973) and vis¬ 
cous (Lightman and Eardley, 1974). A thorough investigation of these was carried 
out by Shakura and Sunyaev (1976) who found that the instabilities develop via 
a bifurcation of unstable modes. Both types of instability occur when radiation 
provides a substantial contribution to the total pressure. In thin discs, the ther¬ 
mal time scales are generally much shorter than the viscous ones and so thermal 
instabilities, if present, dominate over viscous instabilities. The hrst discussions 
of this were made in terms of local linear analysis, in which the global response 
of the disc and nonlinear effects were not taken into account, but then Lightman 
(1974) carried out time-dependent calculations to examine the global behaviour. De¬ 
tailed time-dependent disc evolution was studied hrst in the context of dwarf novae 
(Meyer and Meyer-Hofmeister, 1982; Smak, 1982; Papaloizou, Faulkner and Lin, 
1983). These authors showed that the dwarf nova outbursts arise from collective 
relaxation oscillations of accretion discs in which cool hydrogen ionization regions 
provide a triggering mechanism. As a guide for understanding the global instabili¬ 
ties, an associated local analysis was carried out and a fundamental role in this was 
played by the fact that when the values of the accretion rate, M, and the surface 
density, S, at a particular location (value of r) are plotted against each other for 
a sequence of equilibrium models, the points form characteristic S-shaped curves in 
the log M — log S plane. Analogous investigations of the time-dependent evolution 
of accretion discs around neutron stars and black holes were started by Taam and 
Lin (1984) who found that if the viscous stress tensor is proportional to the total 
pressure, global analysis conhrms the results inferred from local analysis and the in¬ 
stability manifests itself in terms of short time-scale luminosity fluctuations. Similar 
studies were performed subsequently by Lasota and Pelat (1991) and more recently 
by Chen and Taam (1994). Lasota and Pelat concluded that the (locally) thermally 
unstable region in the inner part of the disc always becomes effectively optically 
thin and that the evolution and disc structure depends sensitively on the treatment 
of the transitions between the optically thick and optically thin regimes. Chen and 
Taam (1994) investigated the global evolution of thermal-viscous instabilities in the 
case where the viscous stress is proportional to the gas pressure with a particular 
form for the viscosity parameter. Their calculations show that the instabilities are 
globally coherent throughout the entire unstable region of the disc. 

All of the authors mentioned above studied Keplerian discs. The only investigation 
which we know of in which a full set of time-dependent equations has been solved, 
is the one by Houma, Matsumoto and Kato (1991). 

The importance for accretion discs around black holes of the transonic nature of 
the flow, deviations from Keplerian rotation and non-local cooling by advection, 
was hrst pointed out by Paczyhski and Bisnovatyi-Kogan (1981) who constructed 
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a simple global equilibrium model of a geometrically thin steady disc. A full treat¬ 
ment, including all of these effects, was formulated by Loska (1982) and Muchotrzeb 
and Paczynski (1982). Following their approach, Abramowicz, Czerny, Lasota and 
Szuszkiewicz (1988) calculated a sequence of slim accretion disc models for a wide 
range of accretion rates and found S-shaped curves in logM — logS plane, a result 
which was later conhrmed by Abramowicz, Kato and Matsumoto (1989). This hnd- 
ing recalls the situation for dwarf novae (although the physical processes leading to 
this shape are quite different in the two cases) and this encouraged Honma, Mat¬ 
sumoto and Kato (1991) to study the nonlinear evolution of the thermally unstable 
models to see whether similar behaviour occurred. They found that, at least for 
high viscosities [a = 0.1), a limit-cycle behaviour does indeed occur which is very 
reminiscent of that for dwarf novae. 

The aim of our present study was to examine the global behaviour of the ther¬ 
mally unstable slim disc models in their original formulation (Abramowicz, Czerny, 
Lasota and Szuszkiewicz, 1988) with low viscosity (a = 0.001). The effects of ad- 
vective cooling are included and our results represent a suitable standard reference 
for making comparison with results from forthcoming calculations in which further 
non-local processes will be included which were not present in the original version 
of the slim disc model. 

In Section 2, we discuss the set of basic equations used for our calculations, pointing 
out which assumptions are made and describing the averaging procedure employed. 
The numerical treatment of the equations is presented in Section 3. In Section 4 
the global behaviour of the models is shown and in Section 5 the same models are 
viewed locally. Section 6 contains comments and conclusions. 


2 Basic equations 


In this section, we present the equations used for calculating non-stationary be¬ 
haviour of the original slim disc model and highlight the assumptions involved in 
this treatment. 


2.1 Evolution equations 

We consider here an axisymmetric non-self-gravitating fluid disc formed in the grav¬ 
itational potential <F of a black hole with mass M and use cylindrical polar coor¬ 
dinates (r. If, z) centred on the black hole. The symmetry axis, ; 2 , coincides with 
the rotation axis of the disc. The behaviour of the fluid is governed by the basic 
hydrodynamical equations for conservation of mass, energy and momentum (radial, 
angular and vertical). We study a one-dimensional formulation in which most of the 
equations are written in a vertically integrated form with the vertical structure being 
approximated under the assumptions that hydrostatic balance always holds in the 
vertical direction and that flow properties above and below the equator {z = 0) are 
the same. In this picture, the velocity field is approximated by Vz = 0, Vr = Vr{r,t) 
and = U(p(r,t). Newtonian mechanics is used and general relativistic effects, 
which become signihcant near to the inner edge of the disc, are simulated by using 
the pseudo-Newtonian potential introduced by Paczynski and Wiita (1980): 
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GM 


( 1 ) 


$ = ^_ 

Vr^ + 2;2 — Tq 


Here is the Schwarzschild radius defined by = 2GM/c^, where G and c are 
the gravitational constant and velocity of light respectively. 

Within this framework, the evolution equations can be written in the following form: 
- Mass conservation equation: 




r dr 


{rvr) 


where D/Dt is the Lagrangian derivative given by 


( 2 ) 


D _ d d 
Dt dt~^ ^dr' 

S = Ti{r,t) is the surface density obtained by vertically integrating the volnme 
density p 

rH 

T, = 2 p dz (3) 

Jo 

with H being the half-thickness of the disc, and Vr = Dr/Dt (which is negative for 
an inflow). For the stationary case, the mass conservation eqnation rednces to the 
statement that the local accretion rate is the same at all locations in the disc 

m(r) = —2TirT,Vr = M, (4) 


where m{r) is the total mass internal to the cylindrical radius r, m = dm/dt and 
M is the rate of increase of mass of the black hole. 


- Radial equation of motion: 


DVr 


1 dp 
p dr 



( 5 ) 


where p is the pressure, I = l{r,t) = rv^{r, f) is the specific angular momentum, is 

the value of I for Keplerian motion with = [—r = [GMr /(r — 

Frr and are the viscous forces: 


F = —— (rr ) 

± rr y < rr) 

rp or 


where 



Trr 2 P 


^ dVr _ 
dr r 
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2 2Vr dVr 
W gA'' ^ 

and fi is the shear viscosity coefficient (bnlk and radiative viscosities are not con¬ 
sidered anywhere in the present treatment). The importance of the viscous forces 
Frr and has been pointed out (e.g. by Papaloizou and Stanley, 1986) and they 
have been considered in the more recent studies but we do not include them here 
since they were not part of the original slim disc model. Also, in accordance with 
the original slim disc model, we do not vertically integrate the radial equation of 
motion but, instead, solve it in the equatorial plane. We comment further on this 
in Section 2.2. 


The stationary version of Eq. (5) is 


Vr 


dVr 

dr 


^ ~ ^l) 

p dr r^ 


( 6 ) 


- Azimuthal equation of motion, for the specific angular momentum, 1: 


Dl 

ITt 


rS dr 



rH 


T^rdz 


(7) 


where the rate of shearing is 

dn 

V = (8) 

which, following Shakura and Sunyaev (1973), is often approximated by 

T^r = -OCP (9) 

where a is a (constant) viscosity parameter. This form, which can be derived from 
Eq. (8) in the case of a (stationary) Newtonian Keplerian disc and with the aid of 
several additional approximations, is the one used in the present paper in accordance 
with the original slim disc model. However, a major objective of our future work 
concerns investigation of the extent to which inclusion of a more physical viscosity 
prescription causes differences in non-stationary behaviour with respect to that of 
the “standard” models calculated here. With expression (9), Eq. (7) can be rewritten 
as 


£L = (r^p 

Dt rT, dr ^ 


where P is the vertically-integrated pressure 


( 10 ) 


rH 

P = 2 pdz (11) 

Jo 

For a stationary model, Eq. (10) can be integrated explicitly and, using the boundary 
condition that there is no viscous torque at the black hole horizon, this gives 

Mil - lo) = 2naPr‘^ (12) 
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where Iq is the specihc angular momentum at the inner edge of the disc. 
- Vertical equation of motion: 


Dvz I dp (9$ ^ 

-i =-h- \- F, 

Dt p dz dz 

where the the viscous force Fzz is given by 


Fz 


pdz 


( 13 ) 


with 




2fi d 
3r dr 


{rvr). 


However, F^z is not included in the present treatment. Also, since we are considering 
hydrostatic equilibrium in the vertical direction, the acceleration term vanishes. We 
note that neglecting the vertical acceleration may be a rather poor approximation 
in some parts of the disc, as suggested by two-dimensional studies (Papaloizou and 
Szuszkiewicz, 1994), and investigation of the effect of not neglecting it is one of the 
important issues to be addressed in our forthcoming study. The approach used here 
is exactly the same as for the stationary slim disc models and is discussed further 
in Section 2.2. 


- Conservation of thermal energy: 


~ Qms A Qrad) (14) 

where S is the entropy per unit mass, T is the temperature, Q^s is the rate at 
which heat is generated by viscous friction (including the effects of stresses in both 
the azimuthal and radial directions) 


Qv 



Vr dVr 
r dr 





(15) 


and Qrad is the rate at which heat is lost or gained by means of radiative energy 
transfer (expressed as the sum of terms for the radial and vertical directions) 


Qrad 


1 d 

dT 

d 

dT' 

r dr 


^ d^ 



(16) 


where = AacT^/Snp is the radiative conductivity, a is the radiation constant and 
K is the opacity. In the present treatment, the hrst terms on the right hand side of 
Eqs. (15) and (16) are omitted. Thermal conductivity is not considered here. Using 
the strategy discussed in Section 2.2, the vertically-integrated form of Eq. (14) can 
be written in the following way 

DT _ arT{dn/dr) TF- (4-3/3) T Dp 

'W ~ 0.67(12 - 10.5/3) ~ 0.67phf(12 - 10.5/3) ^ (12 - 10.5/3) 'p^ ^ ^ 
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where p, p and T are the values of these quantities on the equatorial plane, (3 is the 
ratio of the gas pressure to the total pressure and F~ is the radiative flux away from 
the disc in the vertical direction for which we use the expression 


12acT^ 
SnpH ’ 


as discussed in Section 2.2. The stationary form of Eq. (17) is 

M(I+ 67MT^ = 47rrF-. (18) 

dr dr 

The thermodynamic quantities in the equatorial plane are taken to obey the equation 
of state 


p = A;pT + V (19) 

and the opacity is approximated by the Kramers formula for chemical abundances 
corresponding to those of Population I stars 


K = 0.34 X (1 + 6 X lO^V^-^'^). (20) 

Use of this formula represents the only genuine difference between the original slim 
disc formulation and its time-dependent version considered here. The original slim 
disc model used the opacities of Cox and Stewart (1970) stored in the form of a 
table. We have checked that stationary models calculated with these two different 
opacity treatments do not differ signihcantly from each other. 


2.2 Treatment of the vertical structure 


If the disc is not geometrically thick, it is possible to work in terms of vertically- 
integrated variables and to solve an essentially one-dimensional problem. However, 
passing from standard to vertically-integrated variables cannot be done rigorously 
as the vertical structure of the disc is not known in detail. Several approximate 
procedures have been introduced for dealing with this. We describe here the ap¬ 
proach followed in the present treatment and make comparison with the alternative 
approach of Honma, Matsumoto and Kato (1991) (hereafter referred to as HMK). 

Neglecting the 2 :-viscous stress tensor component and the vertical acceleration term, 
Eq. (13) becomes the standard equation of hydrostatic equilibrium: 


dp (9<h 


( 21 ) 
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If one makes the assumption that p and p can be linked (in the vertical direction) 
by the polytropic relation p oc where N is the polytropic index, Eq. (21) can 

be directly integrated (Hoshi, 1977) to give 


= 2{l + N) — , 
po 

(22) 



d(^) =pA^~Jp] ’ 

(23) 


where the subscript 0 denotes equatorial plane quantities and is the Keplerian 

angular velocity given by = [GM/r(r — r^)^] ' . (Hoshi used the Newtonian 
gravitational potential rather than the pseudo-Newtonian one but Eqs. (22) and 
(23) are the same in either case.) HMK used N = 3 for their vertical integrations 
and in this case, = Spo/po but we follow Paczyhski and Bisnovatyi-Kogan 

(1981) and use 


as our equation for H. 




K 



(24) 


In order to calculate the surface density S, the expression (23) can be inserted into 
Eq. (3) and the integral performed directly. For N = 3 one obtains S = {32/33)p^H 
and this was the expression used by HMK. A very similar result is obtained if p 
is taken to be a linear function of with p(0) = po and p{H) = 0. In this case 
S = poH. If, instead, one took the density to be constant in the vertical direction, 
with p{z) = po, then the result would be S = 2poi7 which is rather different. The 
difference between the hrst two expressions is small and, given the approximation 
introduced when using the polytropic law, there is no particular reason to prefer one 
over the other. We use 


S = pqH, (25) 

for conformity with earlier work. A similar discussion applies for the vertically 
integrated pressure P, for which we use 


P = PoH. (26) 

For the radial equation of motion, we use Eq. (5) with the quantities involved tak¬ 
ing their equatorial-plane values. HMK, on the other hand, used the vertically- 
integrated form 


- Ij) P dlnVt^ 

Dt S (9r S dr 


(27) 


which contains the new term {P/lP){dlnVLj^/dr) coming from the vertical integration 
of the gravitational potential. The question of the importance of this term for time- 



dependent calculations has been investigated by Chen (private communication) who 
concludes that it does not make any signihcant difference. 

The thermodynamic relation for the entropy used in our version of the thermal 
energy conservation equation is calculated for a mixture of black body radiation 
and a simple perfect gas with the polytropic index, N, being set equal to 3/2. This 
gives 


T 


DS 


p 

p 


(12 - 10.5/3) 


1 DT 
flDt 


(4-3/3) 


1 Dp 
p Dt 


HMK used the rather similar expression 


(28) 


pT 


DS 

'm 


1 


Dp 


Ta-l 


Dt 


j. P Dp 


(29) 


which is derived for cy = constant (cy being the specihc heat at constant volume). 
Here Ti and Ts are the generalized adiabatic indices. In their integration procedure, 
they used the fact that for iV = 3, /3 does not depend on 2 :. The HMK expression 
for the radiative flux F~ (in the integrated energy equation) is 


_ 16acT^ 

SnpH 

which differs from ours by a numerical factor. Comparison between the stationary 
solutions obtained by Szuszkiewicz (1988) and those obtained by HMK shows that 
the results depend only very weakly (if at all) on the procedure used for the vertical 
integration (Abramowicz, Kato and Matsumoto, 1989). It is not necessarily the case, 
however, that the same will hold for time-dependent calculations. Indeed, Shakura 
and Sunyaev (1976) pointed out that numerical estimates for instability growth 
rates and characteristic wavelengths can be sensitive to the procedure adopted for 
approximating the vertical structure. 


3 The method used for solving the equations numerically 


3.1 The Lagrangian evolutionary code 

The equations presented in Section 2 have been solved numerically, using a La¬ 
grangian hnite difference scheme with standard Lagrangian differencing and grid 
organization. The code was adapted from one developed by Miller and Pantano 
(1990) in a different context. The integration domain extended from a suitable po¬ 
sition in the supersonic part of the flow (for most of the models presented here, this 
was set a.t r ^ 2.5 r^) out to several thousand and was divided into a succession 
of comoving annular zones with each one containing a mass 12% larger than the one 
interior to it. The mass of the innermost zone was determined in accordance with 
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numerical convenience. The mass, m, interior to a given zone boundary was used 
as the radially-comoving independent variable, following the usual practice for La- 
grangian calculations, and the equations of Section 2.1 were rewritten accordingly, 
with radial gradients being transformed using 


A 

dr 


27rrS 


A 

dm 


During the progress of each calculation, a regridding was carried out every time 
that the inner edge of the innermost zone crossed 2.5 in the course of being 
accreted. The innermost zone was then removed from the calculation and all of the 
variables were interpolated onto a new grid having a similar structure (in terms of 
intervals in the independent variable m) to that at the initial time. The interpolation 
(which requires great care in order to avoid introducing destructive inaccuracies) was 
carried out using a piecewise cubic method. With the grid organization which we 
use, adequate spatial resolution is normally obtained with 200 - 300 zones (except 
in the special cases mentioned earlier) and the runs have extended for as many as 
10® - 10® timesteps. The increase in mass of the black hole during the calculations, 
due to accretion of material from the disc, was extremely small (always less than 
one part in 10^® for the models studied here). 

The equations are solved for Vr, r, S, T and I as the main dependent variables. 
(From this point onwards, the thermodynamic parameters p, p and T will always 
be taken to refer to values on the equatorial plane.) Our numerical scheme has 
roughly second-order accuracy in space and in time, the latter being achieved with 
the use of intermediate time-levels in a standard way. S and T are treated as zone- 
centre quantities and computed at the full time-level; r is taken as a zone-boundary 
quantity, computed at the full time-level and Vr and I are taken as zone-boundary 
quantities, computed at the intermediate time-level. The time-step is adjusted in 
accordance with the Courant condition and with two additional constraints on the 
fractional variations of p and T in any single time-step. 

The continuity equation (2) is rewritten in a form which is more convenient for the 
numerical calculations: 


— 0, (30) 

where i is the index of zone boundaries, with i — 1/2 referring to a mid-zone, and 
^i-il 2 is the zone area given by 


Ai-i/2 = 7r(r/ -r/_i). 

The energy equation (17) is solved implicitly for T by means of an iterative proce¬ 
dure, using the secant method. Extrapolation of previous time-step values is used to 
provide initial estimates, and convergence to machine accuracy is typically achieved 
with 4 or 5 iterations. In the viscous heating term, the gradient of the angular 
velocity is substituted by the gradient of the Keplerian angular velocity, an approx¬ 
imation which had been used previously in the stationary calculations in order to 
avoid a serious numerical instability which otherwise occurred there (Muchotrzeb 
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and Paczynski, 1982). In the present time-dependent calculations, we followed the 
same line in order to remain as close as possible to the original slim disc approach. 

When values are known for r, S and T, p, p and H can then be found by solving 
simultaneously the integrated equation of vertical hydrostatic equilibrium (24), the 
solution of the integral for S (25) and the equation of state (19). These equations are 
manipulated algebraically to yield a quadratic equation for p which is then solved 
analytically. The solution for p, p and H is embedded in the iteration loop for T. We 
take this occasion to emphasize that it is necessary to maintain very high levels of 
accuracy throughout the various stages of these calculations as there is a very close 
balance of terms in both the energy equation and the radial equation of motion. 

A simple form of numerical diffusion is introduced when solving the radial equation 
of motion (5) for Vr- In the hnite difference representation of this equation, (ur)” 
((n^)* evaluated at the old time level t"") is substituted by 




(31) 


with the value k = 10“^ being selected on the basis of numerical experiments. The 
reason for introducing this modihcation of the standard scheme will be explained in 
Section 6. 

The inner edge of the grid was placed well within the supersonic part of the flow, 
so that the inner boundary conditions could not affect the evolution external to 
the sonic point. For sub-Eddington stationary models, the sonic radius is near to 
that of the marginally stable orbit at r = Sr^, and so r = 2.5is suitable as an 
initial location for the inner edge of the grid. During the subsequent evolution, a 
check is kept to make sure that the sonic point never gets too close to the grid 
boundary. At every time step, it is necessary to set boundary conditions for Vr 
and I at both the inner and outer edges of the grid. At the inner edge, I is set 
equal to its value at the next grid point, which is a very good approximation since 
the specihc angular momentum is very nearly constant in this region, while Vr is 
calculated from the standard equation (5) with the pressure gradient set to zero 
(which is also an excellent approximation). The outer edge of the grid is set at 
several thousand where there is essentially no change in the variables during the 
time of the calculation and so I and Vr are kept constant there. 


3.2 Initial conditions and perturbations 

As initial conditions for the evolutionary calculations, we used the stationary tran¬ 
sonic disc models constructed by Szuszkiewicz (1988) and described by Abramowicz, 
Czerny, Lasota and Szuszkiewicz (1988). The parameters characterizing these mod¬ 
els are: the mass of the black hole M, the accretion rate M which we measure in 
units of the critical value Me = QAttGM/ cn (the accretion rate for which the lumi¬ 
nosity is equal to T^), and the viscosity parameter, a. For models with accretion 

rates up to the critical one, M/Me is equal to the luminosity of the disc measured 
in units of the Eddington luminosity, = AnGMmpc/aT ~ 10^^{M/Mq), erg/s 
where rup is the proton mass and aT is the Thomson cross section for electron scat¬ 
tering. We will use M/Me and interchangeably for the sub-Eddington and 

Eddington models which are the main focus of our present investigation. All of the 
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models discussed here have the same M and a {M = IOMq and a = 10 but 
M/Me varies over a wide range of values. 

For starting the time-dependent calculations, the initial data from the appropriate 
stationary model is first read in from a hie and transferred onto the hnite-difference 
grid. This is a very delicate procedure since, as mentioned earlier, there is a close 
balance between large terms both in the radial equation of motion and in the en¬ 
ergy equation. The former is particularly sensitive. It is important to transfer the 
minimum possible number of variables onto the grid and then to calculate other 
quantities from the representations of these basic ones on the grid. If this is not 
done, inconsistencies arise which destroy the calculation. Despite the fact that the 
evolutionary code uses hve dependent variables (r, S, T and /), it is only nec¬ 
essary to transfer the hrst four of these from the stationary model since / can then 
be calculated from Eq. (5) on the grid bearing in mind that, for a stationary model, 
Dvr/Dt = Vr dvr/dr. Proceeding in this way turned out to be crucial for the success 
of the calculation. 

The process of transferring the initial data onto the hnite difference grid together 
with the action of numerical noise which inevitably enters during the calculation, 
turn out to be sufficient for introducing suitable generalised perturbations into the 
model to trigger growth of unstable modes which may be present. A lower limit 
on the wavelength of perturbations is, of course, set by the grid spacing but this is 
amply hne enough to include the wavelengths of interest. 


4 Global behaviour of the models 


We have found three different types of global behaviour for the present initially 
stationary slim disc models depending on the accretion rate or, equivalently, on the 
luminosity. Models with L < 0.08 are stable. Models with L between 0.09 and 
1 Ljj develop a violently unstable shock-like feature near to the sonic point, which 
is thought to lead to disruption of the disc, and our calculation then terminates. 
Finally, models with super-Eddington luminosities show progressive slow evolution. 
We concentrate here on the sub-Eddington models, leaving discussion of the super- 
Eddington ones for a forthcoming paper. 

As background for our discussion of the results of the time-dependent calculations, 
it is useful to review briefly the predictions of local stability analysis carried out for 
stationary models. The simplest local stability criterion states that if, at a certain 
radius, the ratio of gas pressure to total pressure, (3, is smaller than 0.4, then the 
disc will be unstable at that radius. We have constructed a sequence of stationary 
models, with a range of luminosities, and key parameters are listed in Table 1. 
The local stability criterion has been applied at every radius for each model and 
the locations of predicted unstable regions have been determined. If, instead, we 
use the rehned criterion derived by Pringle (1976), the predicted unstable regions 
remain essentially the same for all of the models except the one with L = 0.07 
for which the range becomes 6.2 — 9.0r^,. This small difference is due to Pringle’s 
result that when electron scattering does not provide the dominant opacity, the disc 
is thermally stable even when the dominant contribution to the pressure comes from 
radiation. The models with luminosity less than 0.07 are predicted to be stable 
on the basis of local analysis. For L > 0.07L^, the disc can be divided into three 
regions: an inner locally stable zone extending from the sonic point (at rsonic) out 
to the inner edge of the unstable region, the unstable zone (the location of which is 
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given in Table 1) and finally an outer stable zone. 

Following Shakura and Sunyaev (1976), it is possible to give a more extended de¬ 
scription of the stability properties of our stationary models, still on the basis of 
local linear analysis. While this description is somewhat incomplete and has limita¬ 
tions, it is nevertheless quite instructive. For each model, we evaluate the bifurcation 
wavelength, which we will call Amin- Perturbations with wavelengths. A, satisfying 
2H < A < Amin give concentric waves travelling across the surface of the disc while 
for A > Amin there are two branches of growing modes which take the form of stand¬ 
ing waves. For the first of these branches, the perturbation growth rate decreases 
with increasing wavelength and in the limit A 3> 77 this becomes the instability 
discovered by Lightman and Eardley (1974). On the second branch, the growth rate 
increases with increasing wavelength and growth of these modes is due to the ther¬ 
mal instability found by Pringle, Rees and Pacholczyk (1973). The characteristic 
growth time, given by Shakura and Sunyaev, is 


56 - 45/3 - 3/32 

30(0.4 - 0) 


(32) 


with the thermal time scale tth being given by 


1 



The growth time, r, is different for different radii. In Table 1 we give its minimum 
value for each model, Tmin, and the location in the disc where this minimum is 
attained, r{Tmin)- Each of the stationary models listed was evolved forward in time 
with our time-dependent code, as described earlier, and the final column of the 
table shows the (physical) time for which each of these calculations was continued. 
This represents either the time at which the calculation was discontinued because of 
catastrophic instability growth or the time beyond which it was no longer considered 
to be of interest to continue. 

The aim of our time-dependent calculations was to determine which models are 
globally unstable and to see how the instabilities evolve with time. Eirstly, as a 
test of the code, we calculated the evolution of the model with L = 0.01 which is 
stable according to the linear analysis. Without modification, the original basic code 
showed even this model becoming unstable but we found that the instability could 
be removed by introducing a very small amount of additional numerical diffusion 
into the scheme, following the prescription of Eq. (31). The additional diffusion was 
then retained for all of the models studied. We will discuss this further in Section 6. 

Models with L = 0.07 and 0.08 were not found to develop any global instability, a 
result which could already be explained within the context of the local analysis since 
the minimum wavelength Amin-, necessary for the onset of the thermal instability, is 
larger than the width of the unstable region for these models. 

For L in the range from 0.09 to 1 the evolution is terminated by the formation 
of a narrow velocity spike adjacent to the sonic point which, once initiated, grows 
very rapidly in amplitude and disrupts the solution. As a representative example 
of these cases, the global behaviour of the model with L = O.IL^ is illustrated in 
Figures 1, 2, 3 and 4 which show the radial profiles of the radial velocity divided 
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by the sound speed (i.e. the Mach number), the surface density, the temperature 
and the thickness of the disc, at several time-levels during the evolution. An initial 
instability begins to grow gradually in the region predicted to be unstable by the 
local linear analysis and the influence of this then spreads inwards to the sonic point, 
where it triggers a second, more violent, instability (as shown in Figure 5). This is 
a fundamentally non-local process. For this particular model, the perturbation in 
Vrjcg associated with the initial instability extends outwards to beyond 10 and, 
by the time that the run is ended, it has saturated in amplitude and a small region 
of positive (outward) radial velocity has appeared, peaking at ~ Tr^. There are 
corresponding variations in the other parameters, which are in the sense of increasing 
temperature and disc thickness and decreasing surface density in the region of the 
initial instability. The sonic point moves progressively inwards (despite what might 
be inferred from Figure 1). As will be seen from the local views of the global 
behaviour presented in Section 5, the local accretion rate at around Sr^, increases 
more quickly than that at and this leads to an emptying of the inner part of 
the disc. For L = 0.1 the velocity spike near to the sonic point starts to appear 
almost simultaneously with the saturation in amplitude of the initial perturbations 
in Vrjcs and T and, once initiated, the growth occurs very rapidly (in ~ 7 x 10“^ s, 
i.e. on the dynamical timescale) with a shock-like structure developing at the leading 
edge of the spike. There are corresponding rapid changes in the prohles of the other 
variables. The growth of this feature appears to be catastrophic and, certainly, we 
are not able to continue the evolution further with our present numerical scheme. 
However, we are not able to completely rule out the possibility that there might 
be a stabilization in the non-linear regime which could be followed with a more 
sophisticated numerical treatment. We have checked, a posteriori, that even under 
these most extreme conditions the models remain optically thick, according to the 
criterion (/T;efi;jj)^/^(S/2) > 1, both in the expanded innermost part and in the 
region of the velocity spike. Since we strongly suspect that the velocity spike will be 
suppressed or severely altered as a result of including improvements to the original 
slim disc model, we delay further consideration of its origin and behaviour until 
a subsequent paper where the effects of such improvements will be systematically 
discussed. 

While the case L = 0.1 discussed in detail above, is largely representative of 
the others with L in the range from 0.09 to 1 there is a difference appearing for 
models with slightly higher luminosity which needs to be mentioned here. With the 
case just considered, the run is halted by the growth of the velocity spike almost im¬ 
mediately the saturation occurs in the initial perturbation. However, it can be seen 
(particularly in Fig. 2 for T) that the feature associated with the initial instability 
is then just starting to propagate outwards. In higher luminosity models, such as 
the one with L = 0.2 which we will discuss in the next section, this feature has 
time to propagate outwards signihcantly before the run is ended. It is tempting to 
speculate that this might, perhaps, have been the beginning of a cyclic behaviour if 
the instability near to the sonic point had not intervened. 


5 Local view of the global behaviour - S-curves 
and evolutionary tracks 


The existence of S-curves in the log M — log S plane for sequences of models of 
stationary accretion discs around black holes, provided a motivation to look for 
nonlinear cyclical behaviour in these systems in analogy with the situation for dwarf 
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novae. The local stability of a disc at a given radius, r, is related to the location 
of the model on the S-curve for that radius. The lower and middle branches of the 
S-curve correspond to the standard Shakura/Sunyaev models and are simply loci 
representing the balance between viscous heating and radiative cooling. Discs on 
the lower branch are thermally stable while those on the middle branch are unstable 
and a transition via marginal stability occurs at the lower turning point. The upper 
branch has a positive slope which might indicate stability, but here the situation is 
different from that for dwarf novae. On the upper branch here, the heat generated 
by viscosity is mainly advected with the accretion flow, rather then being mainly 
radiated away as on the lower branches, and standard local analysis is not adequate 
for determining whether these models are stable or not. 

Having made our global calculations as described in the previous section, we can 
then take a local view of them, considering the time variation of parameters at a 
particular hxed location r. Following earlier work on local analysis of stationary 
models, it is convenient to consider plots of logrh against logS but it should be 
noted that the trajectories exit from these plots whenever there is a local outflow, 
m < 0. For the case shown in Fig. 1 (L = 0.1 T^), the local outflow is seen only at 
the very end of the run but for models with higher luminosity, the outflow feature 
can sometimes be seen passing by the monitoring location. In Figure 6, we show 
the trajectory of parameters at r = for a model with L = 0.2(Here and 
in the rest of this section, m is measured in units of the critical accretion rate. 
Me ) The dashed lines mark rapid motion out of (or into) the parameter space, 
corresponding to m changing sign. In this case, m remained negative for 14 seconds 
before becoming positive again. 

This hgure is the first of a sequence of similar ones showing different models and 
monitoring locations and so we take this opportunity to describe the way in which 
all of them should be viewed. Stationary models covering the whole range of lu¬ 
minosities are located along the S-shaped curve shown by the dotted line and we 
plot the phase-space trajectory for the particular model of interest at the selected 
monitoring location together with corresponding trajectories for a set of comparison 
models having initial m = 0.01, 0.07, 0.08, 1 and 10. We recall that for the initial 
stationary models, m(r) = M = const and m measured in units of Me is equal to 
LjLj^ for those models which are not super-Eddington. Identification of the differ¬ 
ent initial models is straightforward as it is sufficient to look at the corresponding 
value of m on the axis. For each model, the open square shows the initial state and 
the hlled square shows the state at the end of the calculation with the line joining 
them indicating the phase space trajectory followed. Models for which it is not 
possible to see any change in location on the scale of this plot have the hlled square 
superimposed on the open one so that only the hlled square is seen. Models with 
initial m = 0.01, 0.07 and 0.08 are globally stable and always remain at essentially 
the same point in the phase space for the given radius. The model with Eddington 
luminosity, m = 1, evolved for only a very short time before the run terminated 
and so it is only for small radii (where the evolutionary time-scale is shorter) that 
any movement can be seen in the phase space. Finally the super-Eddington model, 
which we are not discussing in detail here, shows similar sorts of trajectory at the 
diherent radii. 

Although local analysis for discs around neutron stars and black holes is normally 
discussed in terms of the m — S relation, most authors presenting results from time- 
dependent calculations have shown local views in the log T — log S plane instead 
of the logm — logS one. This is reasonable since T and S are primary quantities 
calculated in the numerical codes, whereas m is a secondary quantity used just for 
interpretation, and also there is the problem mentioned earlier, that negative values 
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of m cannot be represented in a plot of logm against logS. In Figure 7 (which 
focuses on our “standard” model with m = 0.1) we show trajectories in both planes 
for comparison. Four different monitoring locations are used (at 3, 5, 8 and 10 r^) 
and the frames are arranged in two columns, one for log m — log S and the other 
for logT — logS. The general conclusions which can be drawn from considering the 
qualitative differences between the evolutionary tracks for different radii, are that 
different parts of the disc do not react in a synchronized way and that the instability 
behaviour seen has a global character. Our discussion will mainly be related to the 
log m — log S curves. 

For r = 3r^, viewed in terms of logm — logS, the model evolves along the lower 
branch of the S-curve for 178 s but then overshoots the position where the S-curve 
bends and makes an attempt to come back again. However, soon after it turns 
and changes its direction, it starts to follow a quite dramatic path and the run is 
halted soon after this. We know, from the global behaviour, that this is caused by 
the development of the instability at the sonic point. The evolutionary track in the 
logT — logS plane is less dramatic, but it is clear also from this that the evolution 
is not going to follow any closed path as would be necessary for having a limit cycle 
behaviour. The situation at 5 if viewed in isolation, would have seemed promising 
for producing a limit cycle behaviour, especially looking at the log T — log S plane 
where the model follows a well-behaved regular path, but it is important to notice 
that the local accretion rate here grows signihcantly less than that at 3 giving rise 
to an unbalanced situation producing a drop in surface density. The local accretion 
rates at 8 and 10 decrease in time, making the effect even more pronounced. This 
description complements the global view of the behaviour of the surface density 
shown in Figure 3. 

For models with higher accretion rates, the importance of the transonic nature of 
the flow, deviation away from Keplerian rotation and advective cooling becomes pro¬ 
gressively greater. To investigate how these non-local effects influence the evolution 
of models with higher accretion rates, we made calculations for several models with 
luminosities between 0.1 and 1 L^. In Figure 8, we show the evolutionary paths for 
L = 0.2, 0.4 and 0.8 at 3 and Sr^. The evolution for L = 0.2 (which we have 
already discussed in connection with Fig. 6, drawn for 8r^) proceeds in a rather 
similar way to that for 0.1 except that the velocity perturbation is here able to 
propagate outwards signihcantly and the run terminates with unstable behaviour at 
the sonic point after a shorter time (71s in this case). These features might seem 
contradictory at hrst sight but it should be noted that the radial velocities are higher 
for 0.2 and so everything proceeds more rapidly. Moving to models with higher 
luminosities, already for L = 0.4 the situation becomes signihcantly diherent. 
Now, the local accretion rate increases systematically at all of the radii considered 
in Fig. 7. The Mach number perturbation propagates out to around 25 but never 
gives any local outhow. For L = 0.8 the Mach number perturbation does not have 
enough time to propagate outwards before the unstable behaviour at the sonic point 
terminates the run (after 25 s). This is directly connected with the fact that the 
drift time t^r = r/vr (which becomes progressively shorter for models with higher 
accretion rates) is now comparable with the thermal time-scale (Szuszkiewicz, 1990). 
By this stage, it is no longer the case that the thermal instability determines the 
evolution of the disc. 

Our results indicate that the models presented here with luminosities in the range 
0.09 - 1 Tg are globally unstable and that the demarcation line between globally 
stable and globally unstable models lies between 0.08 and 0.09It is interesting 
to compare the trajectories in the logT — logS plane for models on either side of 
the demarcation line. Figure 9(a) shows trajectories for two models above it: 0.09 
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and 0.1 at r = 8 r^. We focus attention on the one with 0.09 (which is nearer 
to the demarcation line) and note the evolution starting along an unwinding spiral 
and then diverging off towards higher temperatures. From the corresponding time 
series for the temperature variations, shown in Figure 9(b), we see the temperature 
starting to oscillate with growing amplitude but then quickly turning to unbounded 
growth. Figure 10(a) shows the trajectory in the logT — logS plane for the model 
with L = 0.08at r = '^r^. The evolution proceeds along the inspiralling path 
from the initial, slightly perturbed, state to the asymptotic one with only very small 
changes in the temperature and surface density. (All of this, and the associated 
behaviour in the logm — log S plane, was hidden inside the hlled squares in Figures 
6 , 7 and 8.) The corresponding time series is shown in Figure 10(b) where one can 
see oscillations being damped to zero amplitude. 


6 Comments and conclusions 


In this section, we would hrst like to comment on issues connected with our intro¬ 
duction of additional numerical diffusion. This was motivated by observing that, 
in our hrst calculations, the supposedly stable model with L = 0.01 did not re¬ 
main unchanging but instead developed an instability in the region 3 — 6r^ after 
10 — 20 s of the evolution. This appeared hrst at around 4.7 — 4.8 and grew 
rapidly, causing termination of the calculation. Figure 11 shows the velocity and 
surface-density prohles for the developed instability. This behaviour was not entirely 
unexpected since it was already known that additional global instabilities may be 
introduced as a direct consequence of literal adoption of the Shakura/Sunyaev a 
viscosity prescription (see, for example, Abramowicz, Papaloizou and Szuszkiewicz, 
1993). However, it was important to be sure that the behaviour seen was not merely 
an artefact of our numerical scheme and so we carried out several experiments to 
check on this including using an implicit method for the radial velocity equation 
and a three-point formula for the pressure gradient (the regridding procedure had 
already been subjected to extensive testing). The instability persisted unchanged in 
each case. Another set of experiments was carried out to test the ehect of varying the 
viscosity parameter a. The instability was always present but appeared at different 
locations and after different times. We also tried different viscosity formulations but 
found that the instability was always present. F. Honma (private communication) 
has performed a calculation for the same model with his code but did not hnd the 
same behaviour. Our suspicion is that the discrepancy arises because of the fact 
that our basic code is very non-diffusive and, bearing in mind that the standard 
slim-disc model does not include several potentially important diffusive terms, we 
decided to try adding a low level of numerical diffusion to the radial equation of 
motion, as discussed earlier. We found that the instability was suppressed even 
with extremely small values of k, suggesting that the additional diffusive terms dis¬ 
cussed in Section 2, but not included in the present formulation, probably play a 
crucial role for avoiding this instability which would otherwise make these models 
unsuitable as representations of real astrophysical discs. Since we wanted to focus 
on the more interesting non-stationary behaviour described in the main part of this 
paper, we included the additional diffusion for all of our calculations to avoid them 
being terminated prematurely. Without doing this, all of our sub-Eddington models 
suffer in the same way as the one with L = 0.01 but no similar behaviour is 
seen for super-Eddington models. We emphasize that the kind of numerical scheme 
which we are using here has been extensively tested for a wide range of physical 
situations and does not normally require the addition of extra numerical diffusion 
to give stable behaviour in circumstances such as these. 
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The aim of this paper has been to study the nonlinear evolution of the thermal 
instability for the original version of the slim disc model and to verify predictions 
concerning the stabilizing effect of advection. Setting aside the instability discussed 
in the previous paragraph, we have found that models predicted to be stable by 
local analysis do remain stable and stationary. The same applies for models which, 
according to the local analysis, would have a potentially unstable region smaller 
than the minimum wavelength for unstable perturbations. In terms of luminosities, 
this means that all models with luminosity less than or equal to 0.08 are globally 
stable. For models with luminosities between 0.09 and 1we see the run ended 
by unstable behaviour near to the sonic point. 

For the models studied here and the particular formulation adopted, it seems very 
unlikely that advection alone can be sufficient to allow for the existence of limit cycle 
behaviour. Non-local effects which are not included here (such as radiative diffusion 
in the radial direction) might well play important role in the global evolution of the 
thermal instability since they introduce parabolic terms which could signihcantly 
affect the delicate radial force-balance. 

Our results are rather different from those reported by the Japanese group (HMK). 
However, it should be emphasized that they used a viscosity parameter two orders of 
magnitude larger than ours and this gives rise to models which are very different from 
the original slim disc formulation which we have been studying here. In particular, 
it means that they are dealing with much larger radial velocities and hence the 
effect of advection is much stronger. We are currently in the process of investigating 
models with larger a. 

In our future work, we plan to follow the programme outlined at the beginning of 
the paper: namely, to use the present study as a standard against which to judge 
the effects of systematically adding further sophistications onto the original slim 
disc formulation. Also, we have been studying super-Eddington models. A global 
analysis is really essential for discussing the stability of these since the local approx¬ 
imation is rather poor in this case. However, some qualitative predictions have been 
made on the basis of local analysis (Abramowicz et ah 1988, Wallinder 1991) and 
it will be interesting to compare them with the calculated global behaviour. 
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Table 1: Results from local stability analysis of the initial models 
and the total physical time for which each model was evolved. 


LIL, 

^ sonic 

[ra] 

unstable region 

[re] 

^min 

[ra] 

[rc] 

^min 

N 

tcalc 

N 

0.01 

2.96 

stable everywhere 





0.07 

2.95 

6.2 - 9.7 

28.6 

7.4 

226 

> 5650 

0.08 

2.95 

5.3 - 12.8 

7.7 

6.6 

56.8 

> 7120 

0.09 

2.95 

4.9 - 15.2 

4.7 

6.3 

33.4 

582 

0.1 

2.94 

4.6 - 17.5 

3.6 

5.9 

24.3 

182 

0.2 

2.94 

3.8 - 37.1 

1.6 

4.8 

8.5 

71 

0.4 

2.93 

3.4 - 70.7 

1.1 

4.1 

5.1 

126 

0.5 

2.93 

3.3 - 86.1 

1.0 

3.9 

4.5 

155 

0.6 

2.92 

3.3 - 100.8 

0.9 

3.7 

4.0 

165 

0.7 

2.92 

3.2 - 115.1 

0.8 

3.5 

3.6 

135 

0.8 

2.91 

3.1 - 129.0 

0.7 

3.4 

3.3 

25 

1.0 

2.82 

2.9 - 155.8 

0.5 

3.1 

2.7 

0.08 
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Figure Captions 


Figure 1: Mach number {vr/cg) plotted as a function of r/r^ for the model with 
L = 0.1 at several different times: t = Os (dotted line), t = 150s (short-dashed 
line), t = 162s (long-dashed line), t = 171s (dot/short-dashed line), t = 178s 
(dot/long-dashed line), and t = 180s (solid line). 

Figure 2: Temperature T (in degrees Kelvin) plotted as a function of r/r^ for the 
model with L = 0.1 at several different times. The line description is the same 
as for Figure 1. 

Figure 3: Surface density S (in g/cm^) plotted as a function of r/r^ for the model 
with L = O.IL^ at several different times. The line description is the same as for 
Figure 1. 

Figure 4: The half-thickness of the disc, H (in cm) plotted as a function of r/r^ 
for the model with L = 0.1 at several different times. The line description is the 
same as for Figure 1. 

Figure 5: Development of the instability near to the sonic point {vr/cg = — 1) 
shown in terms of the Mach number. The model concerned is again the one with 
L = 0.1 The dotted line shows the initial Mach number prohle and the rapid 
growth of the sharp-peaked structure is shown by profiles for three consecutive 
output times: t = 180.2001s, 180.2028 s and 180.2069 s. 

Figure 6: The logm — logS plane for r = 8r^ showing the location of the sequence 
of stationary models (dotted line) and the evolutionary paths (solid and dashed 
lines) for models with m = 0.01, 0.07, 0.08, 0.2, 1 and 10. The initial location for 
each model is marked by an open square and the hnal position is marked by a hlled 
square. See text for a more detailed description. 

Figure 7: The left-hand column shows the logm — logS plane for four different 
radii: (a) Sr^, (b) 5r^, (c) 8r^ and (d) 10r^,. The evolutionary paths are for the 
same models as in Fig. 6 except that m = 0.2 is replaced by m = 0.1. The right-hand 
column (hgures (e) - (h)) shows the corresponding behaviour in the logT — logS 
plane. 

Figure 8 : The logm —log S plane for 3r^ (left-hand column) and 5r^ (right-hand 
column) for three different models: (a) m = 0.2, (b) m = 0.4 and (c) m = 0.8. 

Figure 9: (a) The evolutionary paths in the logT — logS plane for r = 8r^ for 
models with rh = 0.09 (solid line) and 0.1 (dashed line). The S-curve is shown by 
the dotted line, (b) The time variation of the temperature at r = 8 for the model 
with 771 = 0.09. 

Figure 10: (a) The evolutionary path in the logT — logS plane for r = 8r^ for 
the model with m = 0.08. The S-curve is shown by the dotted line, (b) The time 
variation of the temperature at r = for the model with m = 0.08. 

Figure 11: The profiles of (a) Mach number and (b) surface density for the model 
with rfi = 0.01 at three different times: t = Os (dashed line), t = 76s (lower 
amplitude perturbation) and t = 87.7s (higher amplitude perturbation). 
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